function y = func_n1_from_p(p, CONSTS)
%     eta = CONSTS.eta;
%     n2 = func_n2_from_p(p, CONSTS);
%     
%     y = - (eta./n2);

    eps = CONSTS.eps;
    g = CONSTS.g;
    q1 = func_q1_from_p(p, CONSTS);
    
    y = -(eps./(p.*g)).*(p.^2+q1.^2+(g^2/eps)-eps);
    
end
